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Nonstoichiometric ceria(Ce02-(5) is a candidate reaction medium to facilitate two step water split- 
ting cycles and generate hydrogen. Improving upon its thermodynamic suitability through doping 
requires an understanding of its vacancy thermodynamics. Using density functional theory(DFT) 
calculations and a cluster expansion based Monte Carlo simulations, we have studied the high tem- 
perature thermodynamics of intrinsic oxygen vacancies in ceria. The DFT+[/ approach was used 
to get the ground state energies of various vacancy configurations in ceria, which were subsequently 
fit to a cluster expansion Hamiltonian to efficiently model the configurational dependence of energy. 
The effect of lattice vibrations was incorporated through a temperature dependent cluster expan- 
sion. Lattice Monte Carlo simulations using the cluster expansion Hamiltonian were able to detect 
the miscibility gap in the phase diagram of ceria. The inclusion of vibrational and electronic entropy 
effects made the agreement with experiments quantitative. The deviation from an ideal solution 
model was quantified by calculating as a function of nonstoichiometry, a) the solid state entropy 
from Monte Carlo simulations and b) Warren-Cowley short range order parameters of various pair 
clusters. 



I. INTRODUCTION 

Ceria(Ce02) and ceria based fluorite structure oxides 
are among the best performing solid oxide fuel cell elec- 
trolytes and have historically been used as automotive 
exhaust catalysts to reduce noxious gases. Due to its 
structural similarity to urania and thoria, ceria is also of 
interest to the nuclear fuels community to investigate ra- 
diation damage to fuel rods in nuclear reactors. Recently, 
Chueh et al.^ demonstrated a two step solar driven ther- 
mochemical cycle based on ceria to split water and pro- 
duce hydrogen as follows. 

Ce02-^, ^ Ce02-5„ + ^^^Oz (1) 

Ce02-^« + {5h - Sl)R20 ^ Ce02-^, + {Sh - <^l)H2 

(2) 

6l and Sh respectively denote the oxygen nonstoichiom- 
etry at the low temperature (T^) and high temperature 
(Th) steps. The high temperature step of the redox cycle 
involves release of oxygen from ceria, forming oxygen va- 
cancies (henceforth referred to as vacancies) in the lattice. 
At the low temperature step, ceria reacts with water, oxi- 
dizing itself and liberating H2 which can be used to gener- 
ate electricity through a fuel cell. From a purely thermo- 
dynamics perspective, these applications stem from the 
remarkable ability of ceria to display significant oxygen 
nonstoichiometry without changing its crystal structure. 
The oxygen vacancy formation is facilitated by the abil- 
ity of Ce to exist in two oxidation states: Ce^+ and Ce^+. 
A general thermodynamic framework based on DFT+ U 
calculations to assess the suitability of an oxide for two- 
step water splitting cycles was proposed by Meredig and 
Wolverton^. 

Experimental work on undoped ceria has been ac- 
tively pursued for a long time, excellent reviews of which 



have been published by Inaba et al."^ and Mogensen et 
al.l^. Thermodynamic data including nonstoichiometry 
(5 as a function of temperature (T) and partial pres- 
sure of oxygen {P02) enthalpy and entropy of 
reduction reaction has been measured extensively^. On 
the computational side, significant research has been de- 
voted to understanding the electronic structure of non- 
stoichiometric ceria from first principles using density 
functional theory(DFT)^ ^ with the standard local den- 
sity approximation (LD A) and generalized gradient ap- 
proximation(GGA) exchange correlation functionals. It 
is well known that the 4f electrons in Ce02 need to be 
treated as valence states to accurately reproduce the ex- 
perimentally observed properties^. Using the conven- 
tional LDA and GGA exchange correlation functionals 
in ceria leads to the self interaction error and a con- 
sequent failure to reproduce the insulating character of 
defect-free ceria. This necessitates adding a Hubbard 
potential([/) correction for the Ce '4f' electrons to gen- 
erate the experimentally observed band gap^. Hybrid 
functionals yield improvement in the electronic picture 
of ceria, but do not significantly change the energetics 
of vacancy formation^. The effects of transition and 
rare earth metal dopant^^^^^ : both aliovalent^^ and 
isovaleniP-^^ on oxygen vacancy formation have been in- 
vestigated using a supercell approach. The oxygen stor- 
age capacity is correlated to the structural relaxation 
brought about by dopants with smaller ionic radius than 
Ce^+ and the electrostatic effects. Activation energies 
for potential vacancy migration pathways have been com- 
puted from first principles ^^""^^ to understand the mech- 
anisms of defect migration at atomistic scale. With the 
development of density functional perturbation theory, 
lattice dynamical properties. Born effective charges and 
phonon density of states have been calculated and t hese 
are found to agree well with experimental dats^^^. 
However, ab initio calculations are currently in- 
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tractable for system sizes greater than a few hundred 
atoms. To obtain finite temperature bulk properties, 
it is necessary to move beyond the realm of DFT cal- 
culations to statistical thermodynamics of larger system 
sizes( >10^ atoms ). While literature abounds with ex- 
perimental investigations of the properties of ceria at 
high temperatures, the computational work described 
above has been primarily limited to studying the elec- 
tronic structure and defect formation energies of ceria at 
absolute zero without attempting to obtain thermody- 
namic properties at higher temperatures, which is where 
its most interesting applications and properties emerge. 

Our computational study aims to isolate and focus on 
the thermodynamics of intrinsic oxygen vacancies in ceria 
relevant to thermochemical cycling. While vacancies are 
central to the thermochemical cycling process, the ther- 
modynamic driving force is governed by the change in S 
(i.e. AS) between Th and Tl. For repeatability of the re- 
dox cycle, it is necessary to operate in the regime of single 
phase nonstoichiometric ceria. Furthermore, at such high 
values of 5, defect interactions become considerable and 
can negatively impact the entropic driving force for the 
reduction of ceria. Understanding the nature of vacancy 
interactions as a function of temperature and concentra- 
tion in ceria will therefore be instrumental in motivating 
dopants to improve its suitability. Aside from its imme- 
diate relevance to thermochemical fuel production, we 
chose ceria as our model system since its phase diagram 
and vacancy thermodynamics are established experimen- 
tally. This would help test the accuracy of the compu- 
tational thermodynamics approach to studying nonstoi- 
chiometric oxides, and extend it to screen dopants and 
predict properties of doped ceria, for which literature is 
not as extensive. 

Modeling of ceria presents a unique challenge: it re- 
quires the ability to capture electron-localization and as- 
sociated electronic entropy effects. A step in this direc- 
tion was made for the Lia;FeP04 system^^ by showing 
that including electronic entropy via a cluster expansion 
approach yields a phase diagram whose topology is in 
qualitative agreement with experiments. We build upon 
this effort and seek to conclusively demonstrate the ther- 
modynamic validity of such an approach, by verifying 
that the inclusion of both electronic and vibrational en- 
tropies results in excellent quantitative agreement with 
experiments for not only the phase diagram, but other 
thermodynamic quantities as well, such as the entropy of 
reduction and short-range order. 

The paper is organized as follows. The methodology 
of studying phase equilibria from first principles as ap- 
plied to ceria will be discussed in Section [TTj Results and 
discussion (section III) will be broken down into sub- 



II. METHODOLOGY 

The standard first principles approach to computing 
phase equilibria has been detailed previously.l^^H^ We 
present a brief overview here to highlight salient fea- 
tures of this approach relevant to nonstoichiometric ox- 
ides. Phase equilibria studies from first principles inte- 
grate both rigorous first principles calculations over se- 
lected small structures and large scale statistical ensem- 
ble based methods. The partition function, which con- 
tains all the thermodynamic information for a system, 
can be coarse grained over a hierarchy of degrees of free- 
dom as described in Equation [3j 



Z= ^^exp [-/3E(a,i 



(3) 



sections to focus independently on electronic structure 
calculations, cluster expansion and free energy integra- 
tion. Wherever applicable, the intricacies of dealing with 
intrinsic oxygen vacancies in ceria will be emphasized. 



where u^a are respectively the vibrational and config- 
urational states of the system constrained to a lattice 
L. Here, the configurational states include both genuine 
configuration variables (the presence of oxygen vacan- 
cies) and electronic state information (the location of the 
Ce^+ ions. The energy E{a^ v) is obtained by performing 
quantum mechanical calculations for a fixed composition 
in the lattice, followed by ionic relaxation. Vibrational 
frequencies of the system are then calculated for small 
displacements away from the relaxed ground state at 
K. The cluster expansion then parametrizes this infor- 
mation in terms of larger structural units and enables 
estimation of energies of cells and compositions inacces- 
sible through first principles calculations in a fast and 
efficient manner. Finally, the thermodynamic integra- 
tion procedure, with the aid of Monte Carlo simulations 
incorporates the effect of compositional fluctuations and 
temperature on the properties of the system and is used 
to derive other thermodynamic quantities. 



A. First principles calculations 

Electronic structure calculations were performed us- 
ing the Vienna ah initio simulation package (VASP), 
a plane wave pseudopotential based DFT code^. GGA 
exchange correlation functional using the PAW (projec- 
tor augmented wave) method^^ was employed. To cor- 
rect for the strong on-site coulombic interaction of the 
Ce 4f electrons, we adopted the rotationally invariant 
GGA+ U formalism introduced by Dudarev et aP^. The 
Hubbard potential term ' V penalizes partial occupancy 
of the f states and opens up a band gap. The value of 
U is typically set by fitting to experimentally established 
band gaps, or quantities such as lattice constant and bulk 
modulus. For nonstoichiometric ceria, based on previous 
LDA+ U and GGA+ U studie^^of oxygen vacancy for- 
mation energies and electron localization on Ce^+, U=b 
for GGA+ U and U =6 for LDA+ U have been proposed 
as optimal. Spin polarized GGA+ U calculations in this 
work were all performed using U =5. 
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Ceria has a fiuorite structure, with oxygen atoms oc- 
cupying the tetrahedral voids of an FCC lattice of Ce. 
Defect calculations were performed on 2x2x2 supercells 
of ceria (96 atoms), using a 2x2x2 A:-point grid. Elec- 
tronic relaxations were performed until the total energy 
difference was less than 10~^ eV while ionic relaxations 
were carried out until residual forces less than 0.02 eV/A 
were achieved. Formation of vacancies in ceria leads to 
expansion of the lattice resulting from increased coulom- 
bic repulsion between the Ce^+ ions and larger ionic ra- 
dius of Ce^+. We account for this by performing multiple 
constant volume relaxations at distinct volumes for each 
structure. The energy benefit accrued from the volume 
relaxation is significant, and if overlooked, could lead to 
erroneous energies predicted with the cluster expansion 
later on. In all, 36 different configurations of vacancies 
were studied with compositions ranging from Ce02 to 
CeOi.75. 

First principles lattice dynamics calculations were per- 
formed to incorporate vibrational effects on phase sta- 
bility. We used the 'small displacement' finite differ- 
ences method as implemented in VASP 5.2 to compute 
the Hessian matrix for the structures. Displacements of 
0.015 A away from the equilibrium relaxed positions were 
employed. For higher vacancy concentrations, we in- 
cluded structures with both pseudo-randomly dispersed 
and clustered arrangements to span the configurational 
dependence of the vibrational frequencies. The force con- 
stants output by VASP were used to obtain the dynami- 
cal matrix at other q-points and calculate out the phonon 
density of states (DOS) using PHONOP^. Gaussian 
smearing and an 8x8x8 q-point mesh were employed for 
the DOS calculation. The vibrational free energy (Fvib) 
and entropy (Svib) were evaluated under the harmonic 
approximation as 



fitting to the database of ab initio energies. 



Sv^b{T) fdFv,^/N 



N 



dT 



(4) 
(5) 



where N is the total number of atoms in the system, u is 
frequency of phonon mode, g{i') is the phonon density of 
states. 



B. Cluster expansion : configurational degrees of 
freedom 



The cluster expansion (CE) Hamiltonian treats config- 
urational disorder by decomposing the energy of a lattice 
into a basis of clusters ( points, pairs, triplets etc.) of 
lattice sites. Each cluster is a polynomial in occupation 
variables and has an associated effective cluster interac- 
tion (ECI), 'J' in equation [6j The ECIs are obtained by 



(6) 



Vacancies are treated as independent species, so any site 
in the anion sub-lattice can be occupied by an oxygen 
ion or vacancy. Additionally, we explicitly treat charge 
state disorder in the cation sub-lattice resulting from 
Ce^+/Ce^+ ( configurational electronic entropy ). In or- 
der to describe the energetics of this system with two 
interacting sub-lattices, we use a multicomponent multi- 
lattice CE formalisni^^^ that works in the product basis 
of cluster functions defined on each sublattice. Despite 
the presence of 4 distinct species (O, Vac, Ce^+, Ce^+), 
constraints of site and charge balance (2[Ce^+] = [Vac]) 
render the system essentially pseudo-binary. Our cluster 
expansion fit was obtained using the mm aps co de in the 
alloy theoretic automated toolkit (ATAT)^^. 

The knowledge of ECIs provides a computationally in- 
expensive and efficient means to compute the energy of a 
large system on the fly during Monte Carlo simulations, 
circumventing the need for time consuming ab initio cal- 
culations. As such, the CE fit to ab initio energies is in- 
dependent of temperature. By cluster expanding phonon 
free energies in the basis of clusters fit to configurational 
energies, the ECIs, and consequently the MC data can 
be made to include vibrational effects. 



C. Lattice Monte Carlo simulations : 
thermodynamic properties 

The fundamental external variables of interest to the 
thermochemical cycling of ceria are temperature (T) and 
oxygen partial pressure (po2 ) • properties are strongly 
dependent on the oxygen nonstoichiometry, which is 
uniquely set for a given (T, ^02)- Thus, ceria lends 
itself to be conveniently studied by semi grand canon- 
ical Monte Carlo simulations, treating temperature and 
chemical potential as external variables. A semi-grand 
canonical ensemble fixes the total number of lattice sites, 
letting the concentration of individual species fluctuate 
in response to an an applied temperature or chemical 
potential change. Consideration of macroscopic charge- 
neutrality leaves us with one independent chemical po- 
tential /i written as : 



/i (/i 



Vac 



M02-) + 2(/ic'e3+ - MCe4+) 



(7) 



where iiyac^l^o'^-il^Ce^+ ^iid iJice'^+ denote the chemical 
potentials of the individual species which are externally 
imposed, ji can be described as the free energy cost as- 
sociated with swapping a pair of Ce^+ and 0^~ with a 
pair of Ce^+ and Vac. 

The Grand Potential <l>(/i,T) with respect to a ref- 
erence can be obtained by thermodynamic integration 
along a fixed T or fixed ji path 

^(To,/i) = $(To,/io)- r {N{To,fi))dfi (8) 
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Kb-L kB-LQ Jto 

(9) 

{E) is the thermodynamically averaged energy, (N) is 
the thermodynamicany averaged concentration, fiQ and 
To are the reference chemical potential and temperature, 
and P = 1/kBT. The low temperature and high temper- 
ature expansions to obtain the reference points for inte- 
gration are not central to this paper and can be found 
elsewhere^^. Our MC runs were performed on 5184 atom 
cells (larger sizes were attempted and found to not af- 
fect the results) using the memc2 code from ATAT^^. 
Temperature steps of 40 K were used for the MC runs, 
with 2000 equilibrium passes and 1000 averaging passes 
at each (T,po2)- Simultaneous spin flips were used to 
maintain charge balance and were constrained to occur 
within two unit cell distances of each other. 

It should be noted that the constraint of charge bal- 
ance alone is not sufficient to guarantee that the sys- 
tem will never undergo a phase separation into multiple 
spurious ground states that are locally non charge bal- 
anced although the overall simulation cell is. This can 
occur when the cluster expansion is only fitted to charge- 
balanced structures, thus providing little guarantee that 
the extrapolated energy of non-charge-balanced struc- 
tures is physically meaningful. We avoided this prob- 
lem by an iterative procedure. Starting with a cluster 
expansion fitted to charge-balanced structure only, we 
monitored the simulation for evidence of phase separa- 
tion into non-charge-balanced structures. Whenever this 
was observed, the energy of the structure was calculated 
from first principles and included in the cluster expan- 
sion to reduce extrapolation errors from the fit. Since 
ah initio calculations imposing periodic boundary con- 
ditions necessarily enforce charge-neutrality, we used a 
neutralizing background charge to estimate the energy of 
non-charge-balanced structures. This ansatz is justified 
whenever the resulting calculated energies are sufficiently 
high, so that these configurations are very rarely sampled 
in equilibrium.) 
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FIG. 1. Phonon density of states from first principles calcu- 
lated under the harmonic approximation. Oxygen vacancies 
lead to softening of vibrational modes in CeOi.91 (3 vacancies 
in a 2x2x2 supercell) compared to Ce02. 



2.5 eV respectively. The formation energy for an oxygen 
vacancy in a 2x2x2 supercell was calculated to be 3.2 eV, 
agreeing with previous GGA+/7 studies^. We included 
up to 8 vacancies in a 2x2x2 supercell (corresponding 
to CeOi.75) in different configurations. Such defect con- 
centrations, while high, have been shown to exist in ce- 
ria under appropriate (T,po2) are of interest to us. 
The lattice expansion associated with vacancy formation 
was evaluated through multiple constant volume relax- 
ations of each configuration and found to be around 2% 
for CeOi.75 (8 vacancies in a 2x2x2 supercell). 

Phonon DOS computed under the harmonic approxi- 
mation are plotted in Fig. [l] Vacancies are clearly stabi- 
lized by vibrations as evidenced by the softening of modes 
in CeOi.91 (3 vacancies in a 2x2x2 supercell) compared 
to Ce02. Further, for a given concentration, clustered 
vacancies had stiffer phonon modes than vacancies dis- 
persed over the supercell. 



B. Cluster Expansion 



III. RESULTS AND DISCUSSION 

A. First principles 

Using GGA+/7 calculations, the equilibrium lattice 
constant for stoichiometric ceria was found to be 5.48A. 
The formation of an oxygen vacancy in a 2x2x2 super- 
cell is accompanied by the localization of two electrons 
onto the 4f states of two Ce atoms with anti-symmetric 
spins. True ground state convergence in other vacancy 
structures was tested by assigning different starting spin 
configurations to the Ce atoms (up or down spin, and 
their locations) and looking for the lowest energy struc- 
ture. The 2p-5d and 2p-4f energy gaps were 5.3 eV and 



The optimal cluster expansion fit to the calculated ab 
initio energies (CE 1) of the 36 relaxed geometries com- 
prised 13 pairs and 1 triplet clusters apart from the null 
and 2 point clusters. Considering that ceria is an ionic 
compound, it is important to ascertain if the long range 
electrostatic interactions are adequately described by a 
short-range cluster expansion. To ascertain this, we fit 
another cluster expansion (CE 2) to the ab initio en- 
ergies after subtracting out the coulombic energy term 
(computed through an Ewald summation). The latter 
was then cluster expanded in the basis of clusters identi- 
fied by the fit, and added as an energetic correction to the 
ECIs. The results are illustrated by Fig. [2] The pair and 
triplet cluster ECIs of CE 2 show identical decay charac- 
teristics with cluster diameter as CE 1, indicating that 
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electrostatic interactions are captured well by the cluster 
expansion. The cross-validation score^^, which provides 
a measure of the predictive power of a cluster expansion 
fit, was close to 0.003 eV for both CE 1 and CE 2. In 
view of these results, we justify using CE 1 for further 
work, given its higher computational efficiency. 

The effect of temperature on ECIs was incorporated 
by cluster expanding the vibrational free energies in the 
basis of clusters of CE 1. This introduced a temperature 
dependent correction to 7 clusters (determined via cross- 
validation) out of a total of 17. 



C. Monte Carlo Simulations 

The Ce-0 phase diagram in the composition range of 
Ce2 03 - Ce02 has been determined by experiments. The 
low temperature portion of the phase diagram is compli- 
cated by the many vacancy ordered phases of composi- 
tion Cen02n-m7 and is not of primary concern to this 
work. Of interest to us is the high temperature phase 
diagram and the ability of ffist principles calculations to 
predict the miscibility gap and vacancy thermodynamics 
in single phase non stoichiometric ceria. 

A temperature-composition plot along a constant /i 
trajectory starting from K is shown in Fig. [Sj Defect 
free Ce02-, the starting ground state (GS), is stable up 
to 1300° C before oxygen vacancies start to form. With 
the inclusion of vibrational effects, the onset of vacancies 
happens at a much lower temperature of 900° C, all other 
variables being the same. Most ah initio phase diagram 
studies have in past neglected the non-configurational 
contributions to defect formation entropy, which can sub- 
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FIG. 2. ECI vs cluster diameter for a cluster expansion fit 
to the as calculated first principles energies (squares, denoted 
by CE 1). Subtracting the electrostatic energy prior to fitting 
and later adding it back (as an energetic correction) does not 
change the variation of ECIs with cluster diameter (filled tri- 
angles, CE 2). This indicates that the electrostatic interac- 
tions are captured by the cluster expansion and do not need 
an explicit treatment. 



stantially affect phase stability. 

In order to accurately predict the miscibility gap, it 
is necessary to perform simulations from two different 
starting points for each /i - a heating simulation from 
the low temperature ordered ground state, and a cooling 
simulation from the high temperature disordered phase. 
At a given temperature, the discontinuity in concentra- 
tion when the grand potentials from the two runs are 
equated leads to the miscibility gap. The phase diagram 
in the composition range of CeOi.8 to Ce02 obtained us- 
ing this approach is shown in Fig. |4] A miscibility gap 
shows up even in the absence of vibrations, but the tem- 
perature scale is off by nearly a factor of two compared 
to experiments. The solubility limit of vacancies in ce- 
ria is underestimated and the miscibility gap is shown to 
persist up to temperatures as high as 1500°C. The tem- 
perature dependent cluster expansion provides a closer 
agreement : the miscibility gap closes at 800° C (690° C 
in experiments) and single phase ceria is shown to be 
stable at higher oxygen nonstoichiometry at any given 
temperature. 

The cluster expansion technique was principally in- 
tended to model configurational disorder in alloys in 
which interatomic interactions tend to be much simpler 
than in insulators or semiconductors. Studies have how- 
ever expanded its domain of applications to describe the 
energetics of Li intercalation in battery materials and 
model charge state disorder through a localized electronic 
entropy term^^, equilibrium composition profile across in- 
terfaces of doped ceria superlattices^^. This is the first 
study of high temperature phase diagram of ceria from 
first principles. That it captures the thermodynamics 
and detects a miscibility gap in a correlated electron sys- 
tem is in itself a significant result; the quantitative agree- 
ment with experiment upon including the effect of vibra- 
tions and electronic entropy shows even more promise. 

The phase diagram helps identify thermodynamically 
stable phases at equilibrium, but even in the region of sin- 
gle phase nonstoichiometric ceria, significant short range 
order can persist leading to deviation of bulk thermody- 
namic properties from that of an ideal solution. Indeed, 
this is the case and ceria deviates from an ideal solution 
model for 5 values as low as 0.007. It is characterized 
experimentally by studying the dependence of nonsto- 
ichiometry (or a dependent property there of, such as 
electrical conductivity)on the partial pressure of oxygen. 
It can also be studied by extracting the nonstoichiom- 
etry dependence of entropy change associated with va- 
cancy formation. Momentarily disregarding the entropy 
of gaseous oxygen released upon vacancy formation, the 
entropy change for the solid phase can be directly eval- 
uated by converting the grand canonical output of MC 
simulations into canonical quantities. 



^-^I^i{ni)) 



(10) 



We fit a model which includes ideal configurational en- 
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tropy with a polynomial correction term (to account for 
non-ideal behavior) to S{6) evaluated per site. 
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where 6 is the nonstoichiometry, A,B,C,D and E are con- 
stants (for a given T). The solid state entropy ?>soiid{^) 
per site referenced to ?iSoiid{^ = 0) is plotted in Fig. p] 
for T=1480 K. The ideal solution entropy(the term whose 
coefficient is A) peaks at (5 = 0.25, but the actual entropy 
of the system plateaus much sooner. The pronounced de- 
viation from ideal solution behavior is apparent from 6 
values as low as 0.01. Vibrations clearly increase the en- 
tropic advantage to having vacancies, but the strongly 
non ideal character persists. To compare with experi- 
mental work^, we computed the entropy of reduction as- 
sociated with forming an oxygen vacancy in the limit of 
infinitesimal change in nonstoichiometry. 
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The solid state entropy is readily available from MC sim- 
ulations. However, there are difficulties associated with 
calculating entropy of gaseous oxygen from first princi- 
ples, mainly concerned with the definition of a reference 
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FIG. 3. A constant fi trajectory MC simulation starting 
from the low temperature ground state (Ce02) illustrating 
vibrational effects on the phase diagram and vacancy concen- 
trations. 
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FIG. 4. (Color online) Miscibility gap in ceria calculated from 
Monte Carlo simulations. The pronounced effect of vibra- 
tions is visible from the suppression of the miscibility gap to 
lower temperatures and enhanced vacancy solubility in non- 



stoichiometric ceria. The experimental phase diagram 
been overlaid for comparison. 
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state. This can be resolved by using the standard state 
entropy data for molecular oxygen from thermochemical 
tables^^. P02 can be obtained using Eqn. 
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fio{T) = ^%{T)^kBT\n 



PO2 



1 atm 



(16) 



While the chemical potentials of the respective species 
from MC simulations are guaranteed to yield the right 
equilibrium composition and free energies, they are ar- 
bitrarily displaced from their true values by an additive 
constant. This prevents a straightforward application of 
Equation [16] As a workaround, we used the experimen- 
tally published nonstoichiometry vs ln(po2) datsP to fit 
the chemical potential from MC and obtain the offset. 
^^Totaii^) computed using this approach at 1480 K is 
illustrated in Fig. [6| The strong deviation from ideal 
solution model suggests that the entropic benefit from va- 
cancies is being offset by some kind of defect interactions. 
These defect interactions could arise from stable defect 
complexes or vacancy ordering over short length scales. 
We look to characterize and quantify the short range or- 
der (SRO) by calculating the thermally averaged pair 
correlations for the pair clusters in real space. Formally, 
this idea is embodied in the Warren- Cowley parameters. 
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For a pair {(0,0,0), (l,m,n)} in the /Vac sub-lattice 
with vacancy site fraction 's', is the occupation vari- 
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FIG. 5. Entropy of Ce02-s referenced to Ce02 computed 
from MC simulations at 1480 K. The deviation from ideal 
solution behavior (configurational disorder in Ce^^/Ce^^ and 
0^~/Vac lattice) becomes apparent even at low 6. Vibrations 
provide entropic benefit to forming vacancies, but the marked 
non-ideal behavior prevails. 
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FIG. 6. Entropy change associated with forming an oxygen 
vacancy (see Eqn. 12) as a function of nonstoichiometry com- 
pared with experiment^ 



able for site i averaged over all sites equivalent by sym- 
metry, (aooo)^ = (1 ~ 2s) ^ is the MC average of the point 
cluster correlation and (<Tooo^/mn) is the averaged correla- 
tion function of the pair cluster. P(cr,T) is the probability 
of configuration a, given by 



P(a,T) = 



Z{{^^i},T) 



exp 



(19) 

where Z({/i^},T) is the semi-grand canonical partition 



fr=aT 




72(110] 



72(111] 




72(100] 



FIG. 7. (Color online) Different vacancy pairs in the anion 
sub-lattice of a unit cell of ceria. Black boxes denote va- 
cancies. The color scheme for the oxygen atoms corresponds 
with that used for plotting SRO parameters of various va- 
cancy pairs in Fig. 8. 



function and {jii} is the set of chemical potentials. If 
vacancies behave ideally (non interacting, randomly dis- 
tributed), then ((Jooo<3-^mn) = (<3-ooo) {crimn) = (1 - 25)^ 
and aimn{^) = for any 6. Clustering of like species is 
given by aimn > and likewise ordering of unlike species 
is given by aimn < 0. 



TABLE I. Cluster type, coordinates and diameter of the first 
five pair clusters in the anion sublattice of ceria. There are 2 
distinct 1/2 [111] clusters, one within the unit cell(superscript 
'a') and another extending out of the unit cell, with a face 
centered Ce atom halfway in between (superscript 'b'). 



Cluster type 


Sites (fractional 
coordinates) 


Diameter(A) 


1/2[100] 


(0.25,0.25,0.25) 
(0.25,0.25,0.75) 


2.74 


1/2[110] 


(0.25,0.25,0.25) 
(-0.25,0.25,0.75) 


3.87 


1/2[111]" 


(0.75,0.75,0.75) 
(0.25,0.25,1.25) 


4.75 


1/2[111]' 


(0.75,0.75,0.75) 
(0.25,0.25,0.25) 


4.75 


[100] 


(0.75,0.75,0.75) 
(-0.25,0.75,0.75) 


5.48 



Fig. [7| shows three vacancy pairs in the anion sublat- 
tice of a unit cell of ceria. There are 2 distinct 1/2 [111] 
clusters, one within the unit cell and another extend- 
ing out of the unit cell, with a face centered Ce atom 
halfway in between. Tab. [T| summarizes the coordinates 
and diameters of the five smallest pair clusters in the an- 
ion sublattice. In Fig. 8(a), a{5) at 1320 K for the five 
clusters in Table |l] are plotted. Close to stoichiometry, 
o^imn = for all (Imn), as would be expected for non 
interacting vacancies. However, deviations from zero be- 
come apparent even for slightly off-stoichiometric com- 
positions. In particular, there is strong preference for va- 
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cancies to align along 1/2 [110]. The negative correlation 
along 1/2 [100] indicates that two vacancies are not ther- 
mo dynamically favored at nearest neighbor sites. The 
two distinct 1/2 [111] pairs have nearly the same value of 
a > and hence are possible directions for short range 
clustering of vacancies. The effect of temperature is to 
favor disorder (Fig. 8(b) for WC parameters at 1791 K), 
as can be seen from the fact that aimn ~ up to larger 
nonstoichiometry. At still higher concentrations of va- 
cancies, the aimn start to deviate from zero and show 
similar clustering/ordering tendencies along the respec- 
tive directions. 
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FIG. 8. (Color online) Warren-Cowley short range order pa- 
rameters (aimn) as a function of nonstoichiometry for various 
pair clusters in the 0^~/Vac sublattice at (a) 1262 K and (b) 
1791 K. Of note is a strong clustering tendency along 1/2 [110]. 



There are multiple related system properties associated 
with vacancy formation in ceria, hence no one rule can 
govern the choice of dopants. While vacancies are crucial, 



it is the change in nonstoichiometry between the high and 
low temperature steps (Eqn. [T]) that ultimately estab- 
lishes the oxygen uptake and consequently the amount 
of hydrogen produced. An important implication of this 
study is that interactions between vacancies (can include 
dopant-vacancy interactions in doped ceria) markedly di- 
minish the entropic driving force for the thermochemi- 
cal cycling of ceria. A dopant which tends to increase 
the dielectric constant of ceria would screen vacancies 
from seeing each other and inhibit short range ordering. 
A computational study combining accurate first princi- 
ples calculations of electronic structure and energetics, 
together with an efficient cluster expansion based Monte 
Carlo simulations of configurational disorder, tempera- 
ture and oxygen chemical potential effects is essential to 
screen potential dopants for ceria. 



IV. CONCLUSIONS 

We successfully employed a cluster expansion Hamilto- 
nian based lattice Monte Carlo simulations approach to 
quantitatively compute the high temperature thermody- 
namics of oxygen vacancies in ceria from first principles. 
The ground state energy and electronic structure of non- 
stoichiometric ceria were obtained from GGA+ U super- 
cell calculations. The database of structures and energies 
was used to fit a coupled cluster expansion that explicitly 
accounts for charge state disorder (Ce^+/Ce^+). Lattice 
vibrational free energies were calculated from first prin- 
ciples under the harmonic approximation and found to 
favor formation of vacancies. Vibrational effects were in- 
corporated as a temperature correction to the ECIs. The 
phase diagram obtained from lattice Monte Carlo simu- 
lations was found to exhibit a miscibility gap. The inclu- 
sion of vibrations resulted in quantitative corrections to 
the composition and temperature range of the miscibil- 
ity gap, yielding excellent agreement with experiments. 
The solid state entropy change resulting from vacancy 
formation was evaluated and the deviation from ideal so- 
lution behavior illustrated through composition depen- 
dence of entropy. To further quantify the defect inter- 
actions leading to deviations from ideality. Warren Cow- 
ley short range order parameters were computed. It was 
found that there is a strong preference for vacancies to 
cluster along 1/2 [110] and 1/2 [111] directions, while the 
nearest neighbor 1/2 [100] sites exhibited ordering behav- 
ior. While temperature does disorder the structure, the 
aforementioned behavior was shown to persist at temper- 
atures as high as 1780 K. 
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